Assuming we have some extra time, this additional exercise is going to explore the species location data in more detail. We will produce two interactive map examples that you can open in your web browser. These types of simple web maps and applications allow you to share your data quickly and easily with project partners, colleagues, and the general public, with a lot of control over style and content. These examples are also meant to emphasize the amazing capabilities of existing APIs that you can leverage to access high resolution imagery and data visualization functionality that would otherwise be cost prohibitive to purchase or code from scratch.
We’ve already loaded a few of these packages in the previous example, but let’s include everything we need for this example just to be safe.
# Packages for Leaflet Mapping Example
library(rgdal)
library(raster)
library(sp)
library(sf)
library(dplyr)
library(tidyr)
library(ggplot2)
library(leaflet)
library(leaflet.extras)
library(leaflet.opacity)
library(htmltools)
library(htmlwidgets)
# Packages for 3D Mapping Example
library(geojsonio)
library(mapdeck)
For this example, we will re-import all of the data and keep all of the species points to make our interactive map.The next few cells will be similar to what you already completed in the earlier exercise.
# Set the working/base directory where our project files are stored. This will be unique to your computer.
baseDir <- getwd()
# Create an output directory for any processed files
dir.create(path = "output")
outDir <- paste0(baseDir,"/output")
Again, we’ve done this, so we will do it in a single cell.
# 1. Read in the species sample points again...
d1 <- read.csv(paste0(baseDir,"/data/cucurbitaData.csv"))
# 2. Read in the protected areas raster again...
protLands <- raster::raster(x = paste0(baseDir,"/data/proAreas.tif"))
# 3. Convert the .csv to a spatial points data frame, again...
coordinates <- d1[,c(4,3)]
proj4 <- protLands@crs
sp1 <- sp::SpatialPointsDataFrame(coords = coordinates,
data = d1,
proj4string = proj4)
# 4. Crop raster to spatial data frame extent
protCrop <- raster::crop(x = protLands, y = sp1)
Next we are going to extract information from the protected areas raster to a new column in our data frame called “protected”. It will fill the each row with a “1” if it falls within a protected area, and an “NA” if it falls outside of one.
sp1$protected <- raster::extract(protCrop, sp1)
head(sp1)
## X taxon latitude longitude type protected
## 1 1 Cucurbita_cordata 28.95470 -113.5625 G 1
## 2 2 Cucurbita_cordata 24.28577 -111.0110 H NA
## 3 3 Cucurbita_cordata 25.93470 -111.5421 H 1
## 4 4 Cucurbita_cordata 25.94749 -111.7300 H 1
## 5 5 Cucurbita_cordata 25.95765 -111.4609 H 1
## 6 6 Cucurbita_cordata 26.20671 -111.6175 H 1
Now we want to label points that fall within protected areas so we can map protected species locations. Since we already extracted protected area information to our data frame, records that don’t fall within protected areas will have an “NA” value. In this example we are going to change the “NA” values to “0”. Another option would be to simply omit all “NA” values from the table if we only wanted the protected points, but in this case we want to visualize the differences.
sp1.class <- sp1
indexNA <- which(is.na(sp1.class$protected))
sp1.class[indexNA, "protected"] <- 0
head(sp1.class)
## X taxon latitude longitude type protected
## 1 1 Cucurbita_cordata 28.95470 -113.5625 G 1
## 2 2 Cucurbita_cordata 24.28577 -111.0110 H 0
## 3 3 Cucurbita_cordata 25.93470 -111.5421 H 1
## 4 4 Cucurbita_cordata 25.94749 -111.7300 H 1
## 5 5 Cucurbita_cordata 25.95765 -111.4609 H 1
## 6 6 Cucurbita_cordata 26.20671 -111.6175 H 1
# We could also omit ALL "NA" values for an "all protected" data frame. But we can ignore this code for now.
# sp1.nona@data <- na.omit(sp1.nona@data)
# head(sp1.nona)
Next we are going to load the species range maps, or species distribution model rasters found in your “/data” folder. We didn’t use these files in the previous example, but they are useful for mapping species ranges along with the sample point data. The file is in “.rda” format, so we will use the load() function and then pull each raster out of the file as individual raster objects.
# Load the .rda file
path <- paste0(baseDir,"/unzip/CucurbitaRasters.rda")
load(path)
# Extract each species' raster to individual objects
cordata <- CucurbitaRasters$cordata
palmata <- CucurbitaRasters$palmata
digitata <- CucurbitaRasters$digitata
# Plot an example...
plot(digitata)
We can see from the map above that the data contains presence and absence values for the range map. We only want to plot the presence data (i.e., “1”) in our web map, so we need to filter the data.
digi.filter <- digitata
digi.filter[digi.filter[] < 1] = NA
digi.filter
## class : RasterLayer
## dimensions : 144, 150, 21600 (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667 (x, y)
## extent : -124.5, -99.5, 20.33333, 44.33333 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## source : memory
## names : digitata
## values : 1, 1 (min, max)
# Check the output...
plot(digi.filter)
Let’s do that for the other species range maps…
corda.filter <- cordata
corda.filter[corda.filter[] < 1] = NA
palma.filter <- palmata
palma.filter[palma.filter[] < 1] = NA
The previous tutorial had us create a simple Leaflet web map as an example, but Leaflet has so much to offer! The following cell will create a fully functional web application, where we can view both point vector and raster layers, make a heatmap of location density, and visualize different base maps. This section is meant to provide a deeper overview of the capabilities of Leaflet for making interactive web maps and applications. Again, one exciting, but very easy feature to add is the ability to easily make heatmaps based on the point location densities.
A great place to choose color values (HEX, RGB, CMYK), especially for finding colorblind safe and print optimized color values, is from is Color Brewer
Here we are grabbing Leaflet provider basemap tiles. You can see all of the options here.
There are many (though not unlimited) options for customizing your web map by utilizing all the features in the Leaflet package. You can extend the map’s functionality with a growing ecosystem of plugins and add-on packages (e.g., leaflet.extras), or just join a future tutorial on JavaScript.
pal <- colorFactor(c("#5ab4ac", "#d95f02", "#7570b3"), domain = c("Cucurbita_cordata", "Cucurbita_palmata", "Cucurbita_digitata"))
propal <- colorFactor(c("#e31a1c", "#1f78b4"), domain = c(0, 1))
map <-
sp1.class %>%
leaflet() %>%
addProviderTiles(providers$CartoDB.Positron, group = "CartoDB", options = providerTileOptions(minZoom = 4, maxZoom = 12)) %>%
addProviderTiles(providers$OpenTopoMap, group = "Topographic", options = providerTileOptions(minZoom = 4, maxZoom = 12, opacity = 0.45)) %>%
addProviderTiles(providers$Esri.WorldImagery, group = "Satellite", options = providerTileOptions(minZoom = 4, maxZoom = 12)) %>%
addCircleMarkers(color = ~pal(taxon), radius = 4, weight = 1, fillOpacity = 0.60, group = "Species Records") %>%
addCircleMarkers(color = ~propal(protected), radius = 4, weight = 1, fillOpacity = 0.60, group = "Protected Status") %>%
addHeatmap(lng=~as.numeric(longitude),
lat=~as.numeric(latitude),
radius = 8,
group = "Density (Heatmap)") %>%
addRasterImage(protCrop, colors = "#fdbf6f", opacity = 0.75, group = "Protected Areas") %>%
addRasterImage(corda.filter, colors = "#1b9e77", opacity = 0.65, group = "Cordata Range") %>%
addRasterImage(digi.filter, colors = "#d95f02", opacity = 0.65, group = "Digitata Range") %>%
addRasterImage(palma.filter, colors = "#7570b3", opacity = 0.65, group = "Palmata Range") %>%
addLegend("bottomleft", pal = pal, values = ~taxon,
title = "Taxonomic Name",
opacity = 1, group = "Species Records") %>%
addLegend("bottomleft", pal = propal, values = ~protected,
title = "Protected Status",
opacity = 1, group = "Protected Status") %>%
addLayersControl(
baseGroups = c("CartoDB (default)", "Topographic", "Satellite"),
overlayGroups = c("Species Records", "Density (Heatmap)", "Protected Status", "Protected Areas", "Cordata Range", "Digitata Range", "Palmata Range"),
options = layersControlOptions(collapsed = TRUE)) %>%
# Set the initial view and zoom level of the map
setView(-115, 30.75, zoom = 5) %>%
# Using hideGroup gives us the option to hid map layers on launch to declutter our map
hideGroup("Density (Heatmap)") %>%
hideGroup("Protected Status") %>%
hideGroup("Protected Areas") %>%
hideGroup("Cordata Range") %>%
hideGroup("Digitata Range") %>%
hideGroup("Palmata Range") %>%
# This button add on allows us to reset the map to the original extent
addResetMapButton()
map
Finally, we can save our map as an .html file so we can explore the map in a web browser!
saveWidget(map, file= paste0(outDir,"leafletMap.html"))
Check your “/outputs” folder for the .html file. Double click and see how the map looks in your web browser of choice.
Okay, one more fun map experiment! Let’s use the mapdeck package to make a 3D web visualization of our data. This package uses “Deck.gl”, a user-friendly WebGL javascript library that integrates nicely with Mapbox, which, like Leaflet, is a great web mapping platform.
WebGL is a JavaScript API for rendering interactive 2D and 3D graphics in your web browser. Make sure your browser is WebGL compatible
WebGL Platforms worth exploring:
Here I want to extract the species sample points to ecoregions so that we have something informative to plot. It’s completely possible to do this with our available data in R, but we only have 15 minutes, so I’ve already prepared a shapefile in ArcGIS Pro with species location attributes added. We haven’t loaded a shapefile into R yet, so let’s do that and quickly view the data.
# Import and view
ecoCount <- rgdal::readOGR(paste0(outDir,"/eco_feature_count_04162021.shp"))
## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\niko8816\Desktop\NK\Jobs\CSU_Geospatial_Centroid\Research-and-Program-Coordinator-main\output\eco_feature_count_04162021.shp", layer: "eco_feature_count_04162021"
## with 18 features
## It has 4 fields
plot(ecoCount,
col=ecoCount$COUNTS,
main = "Species Sample Count By Ecoregion")
# Convert shapefile to GeoJSON
ecoCount.json <- geojson_json(ecoCount)
ecoCount.sf <- geojson_sf(ecoCount.json)
ecoCount.sf
## Simple feature collection with 18 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -123.1593 ymin: 20.41155 xmax: -99.53937 ymax: 42.06112
## Geodetic CRS: WGS 84
## First 10 features:
## ECO_NAME COUNTS latitude longitude
## 1 Apache Highlands 88 31.56698 -111.0514
## 2 Arizona-New Mexico Mountains 2 32.96590 -108.5767
## 3 Baja California Desert 96 24.28577 -111.0110
## 4 California Central Coast 13 35.60691 -120.1635
## 5 California South Coast 340 31.66667 -115.9500
## 6 Chihuahuan Desert 31 32.23387 -108.0502
## 7 Colorado Plateau 1 35.82304 -113.5702
## 8 Great Basin 1 38.10881 -119.0518
## 9 Great Central Valley 62 35.21665 -118.7277
## 10 Gulf Of California Xeric Scrub 38 28.95470 -113.5625
## geometry
## 1 MULTIPOLYGON (((-109.0889 2...
## 2 MULTIPOLYGON (((-105.8239 3...
## 3 MULTIPOLYGON (((-111.6571 2...
## 4 MULTIPOLYGON (((-120.2269 3...
## 5 MULTIPOLYGON (((-119.373 34...
## 6 MULTIPOLYGON (((-99.57135 2...
## 7 MULTIPOLYGON (((-112.3335 3...
## 8 MULTIPOLYGON (((-113.1648 3...
## 9 MULTIPOLYGON (((-118.7917 3...
## 10 MULTIPOLYGON (((-110.3555 2...
Now let’s plot some data to see how it looks. First we will test just using a simple data frame of species sample locations.
This package uses the Mapbox API for the basemap, so we will need to register for a free Mapbox account to access our user token. I’ve included my public token for this example, but please be kind and create your own for any future exploration.
# Set the Mapbox key
key <- set_token("pk.eyJ1Ijoia290bGluaWMiLCJhIjoiZ3FCVFdTZyJ9.WfWEgjDKl9N0NKUyeLlAgA")
mapdeck(
style = mapdeck_style('dark'),
pitch = 45,
location = c(-112, 31),
zoom = 5) %>%
add_grid(
data = d1,
lat = "latitude",
lon = "longitude",
cell_size = 5000,
elevation_scale = 200,
layer_id = "taxon",
update_view = FALSE
) %>%
add_heatmap(
data = d1,
lat = "latitude",
lon = "longitude",
weight = "X",
colour_range = colourvalues::colour_values(1:6, palette = "viridis"),
update_view = FALSE
)
# We are going to exaggerate our data values so that the features "pop" more.
ecoCount.sf$COUNTS <- ecoCount.sf$COUNTS * 100
mapdeck(
style = mapdeck_style('dark'),
pitch = 45,
location = c(-112, 31)
) %>%
add_polygon(
data = ecoCount.sf,
layer = "polygon_layer",
fill_colour = "ECO_NAME",
elevation = "COUNTS",
legend = TRUE
)
ecodeck <- mapdeck(token = key, style = mapdeck_style("dark")) %>%
add_polygon(
data = ecoCount.sf,
layer = "polygon_layer",
fill_colour = "ECO_NAME",
elevation = "COUNTS",
legend = TRUE
)
# We are going to save the widget as an .html file again so we can open it in our browser.
saveWidget(ecodeck, file="3Dmap.html")
Good work! You’ve finished this additional portion of the tutorial.